Introduction à la modélisation statistique bayésienne

Planning

Cours n°01 : Introduction à l’inférence bayésienne
Cours n°02 : Modèle Beta-Binomial
Cours n°03 : Introduction à brms, modèle de régression linéaire
Cours n°04 : Modèle de régression linéaire (suite)
Cours n°05 : Markov Chain Monte Carlo
Cours n°06 : Modèle linéaire généralisé
Cours n°07 : Comparaison de modèles
Cours n°08 : Modèles multi-niveaux (généralisés)
Cours n°09 : Examen final

\[\newcommand\given[1][]{\:#1\vert\:}\]

Langage de la modélisation

\[ \begin{align} y_{i} &\sim \mathrm{Normal}(\mu_{i}, \sigma) \\ \mu_{i}&= \alpha + \beta x_{i} \\ \alpha &\sim \mathrm{Normal}(60, 10) \\ \beta &\sim \mathrm{Normal}(0, 10) \\ \sigma &\sim \mathrm{HalfCauchy}(0, 1) \end{align} \]

Objectif de la séance : comprendre ce type de modèle.

Un premier modèle

library(tidyverse)
library(imsb)

d <- open_data(howell)
str(d)

'data.frame':   544 obs. of  4 variables:
 $ height: num  152 140 137 157 145 ...
 $ weight: num  47.8 36.5 31.9 53 41.3 ...
 $ age   : num  63 63 65 41 51 35 32 27 19 54 ...
 $ male  : int  1 0 0 1 0 1 0 1 0 1 ...

d2 <- d %>% filter(age >= 18)
head(d2)

   height   weight age male
1 151.765 47.82561  63    1
2 139.700 36.48581  63    0
3 136.525 31.86484  65    0
4 156.845 53.04191  41    1
5 145.415 41.27687  51    0
6 163.830 62.99259  35    1

Un premier modèle

\[h_{i} \sim \mathrm{Normal}(\mu, \sigma)\]

d2 %>%
    ggplot(aes(x = height) ) +
    geom_histogram(aes(y = ..density..), bins = 20, col = "white") +
    stat_function(fun = dnorm, args = list(mean(d2$height), sd(d2$height) ), size = 1)

Loi normale

\[ p(x \given \mu, \sigma) = \frac{1}{\sqrt{2 \pi \sigma^{2}}} \exp \bigg[-\frac{1}{2 \sigma^{2}} (\mu - x)^{2} \bigg] \]

data.frame(value = rnorm(n = 1e4, mean = 10, sd = 1) ) %>% # 10.000 samples from Normal(10, 1)
    ggplot(aes(x = value) ) +
    geom_histogram(col = "white")

D’où vient la loi normale ?

Contraintes : Certaines valeurs soient fortement probables (autour de la moyenne \(\mu\)). Plus on s’éloigne, moins les valeurs sont probables (en suivant une décroissance exponentielle).

D’où vient la loi normale ?

\[ y = \exp \big[-x^{2} \big] \]

On étend notre fonction aux valeurs négatives.

D’où vient la loi normale ?

\[ y = \exp \big[-x^{2} \big] \]

Les points d’inflection nous donnent une bonne indication de là où la plupart des valeurs se trouvent (i.e., entre les points d’inflection). Les pics de la dérivée nous montrent les points d’inflection.

D’où vient la loi normale ?

\[ y = \exp \bigg [- \frac{1}{2} x^{2} \bigg] \]

Ensuite on standardise la distribution de manière à ce que les deux points d’inflection se trouvent à \(x = -1\) et \(x = 1\).

D’où vient la loi normale ?

\[ y = \exp \bigg [- \frac{1}{2 \color{steelblue}{\sigma^{2}}} x^{2} \bigg] \]

On insère un paramètre \(\sigma^{2}\) pour contrôler la distance entre les points d’inflection.

D’où vient la loi normale ?

\[ y = \exp \bigg [- \frac{1}{2 \color{steelblue}{\sigma^{2}}} (\color{orangered}{\mu} - x)^{2} \bigg] \]

On insère ensuite un paramètre \(\mu\) afin de pouvoir contrôler la position (la tendance centrale) de la distribution.

D’où vient la loi normale ?

\[ y = \frac{1}{\sqrt{2 \pi \color{steelblue}{\sigma^{2}}}} \exp \bigg[-\frac{1}{2 \color{steelblue}{\sigma^{2}}} (\color{orangered}{\mu} - x)^{2} \bigg] \]

Mais… cette distribution n’intègre pas à 1. On divise donc par une constante de normalisation (la partie gauche), afin d’obtenir une distribution de probabilité.

Modèle gaussien

Nous allons construire un modèle de régression, mais avant d’ajouter un prédicteur, essayons de modéliser la distribution des tailles.

Modèle gaussien

On définit \(p(\mu, \sigma)\), la distribution a priori conjointe de tous les paramètres du modèle. On peut spécifier ces priors indépendamment pour chaque paramètre, sachant que \(p(\mu, \sigma) = p(\mu) p(\sigma)\).

\[\color{steelblue}{\mu \sim \mathrm{Normal}(178, 20)}\]

Modèle gaussien

On définit \(p(\mu, \sigma)\), la distribution a priori conjointe de tous les paramètres du modèle. On peut spécifier ces priors indépendamment pour chaque paramètre, sachant que \(p(\mu, \sigma) = p(\mu) p(\sigma)\).

\[\color{steelblue}{\sigma \sim \mathrm{Uniform}(0, 50)}\]

Visualiser le prior

library(ks)
sample_mu <- rnorm(1e4, 178, 20) # prior on mu
sample_sigma <- runif(1e4, 0, 50) # prior on sigma
prior <- data.frame(cbind(sample_mu, sample_sigma) ) # multivariate prior
H.scv <- Hscv(x = prior, verbose = TRUE)
fhat_prior <- kde(x = prior, H = H.scv, compute.cont = TRUE)
plot(
    fhat_prior, display = "persp", col = "steelblue", border = NA,
    xlab = "\nmu", ylab = "\nsigma", zlab = "\n\np(mu, sigma)",
    shade = 0.8, phi = 30, ticktype = "detailed",
    cex.lab = 1.2, family = "Helvetica")

Prior predictive checking

sample_mu <- rnorm(1000, 178, 20)
sample_sigma <- runif(1000, 0, 50)

data.frame(x = rnorm(1000, sample_mu, sample_sigma) ) %>%
    ggplot(aes(x) ) +
    geom_histogram() +
    labs(x = "Taille (en cm)", y = "Nombre d'échantillons")

Fonction de vraisemblance

mu_exemple <- 151.23
sigma_exemple <- 23.42

d2$height[34] # une observation de taille (pour exemple)

[1] 162.8648

Fonction de vraisemblance

On veut calculer la probabilité d’observer une certaine valeur de taille, sachant certaines valeurs de \(\mu\) et \(\sigma\), c’est à dire :

\[ p(x \given \mu, \sigma) = \frac{1}{\sqrt{2 \pi \sigma^{2}}} \exp \bigg[-\frac{1}{2 \sigma^{2}} (\mu - x)^{2} \bigg] \]

dnorm(d2$height[34], mu_exemple, sigma_exemple)

[1] 0.01505675

Fonction de vraisemblance

\[ p(x \given \mu, \sigma) = \frac{1}{\sqrt{2 \pi \sigma^{2}}} \exp \bigg[-\frac{1}{2 \sigma^{2}} (\mu - x)^{2} \bigg] \]

Ou avec une fonction maison…

normal_likelihood <- function (x, mu, sigma) {
  
  bell <- exp( (- 1 / (2 * sigma^2) ) * (mu - x)^2 )
  norm <- sqrt(2 * pi * sigma^2)
  
  return(bell / norm)
  
}

normal_likelihood(d2$height[34], mu_exemple, sigma_exemple)

[1] 0.01505675

Distribution postérieure

\[ \color{purple}{p(\mu, \sigma \given h)} = \frac{\prod_{i} \color{orangered}{\mathrm{Normal}(h_{i} \given \mu, \sigma)}\color{steelblue}{\mathrm{Normal}(\mu \given 178, 20)\mathrm{Uniform}(\sigma \given 0, 50)}} {\color{green}{\int \int \prod_{i} \mathrm{Normal}(h_{i} \given \mu, \sigma)\mathrm{Normal}(\mu \given 178, 20)\mathrm{Uniform}(\sigma \given 0, 50) \mathrm{d} \mu \mathrm{d} \sigma}} \]

\[ \color{purple}{p(\mu, \sigma \given h)} \propto \prod_{i} \color{orangered}{\mathrm{Normal}(h_{i} \given \mu, \sigma)}\color{steelblue}{\mathrm{Normal}(\mu \given 178, 20)\mathrm{Uniform}(\sigma \given 0, 50)} \]

Distribution postérieure - Grid approximation

# on définit une grille de valeurs possibles pour mu et sigma
mu.list <- seq(from = 140, to = 160, length.out = 200)
sigma.list <- seq(from = 4, to = 9, length.out = 200)

# on étend la grille à toutes les combinaisons possibles de mu et sigma
post <- expand.grid(mu = mu.list, sigma = sigma.list)

# calcul de la log-vraisemblance (pour chaque combinaison de mu et sigma)
post$LL <-
  sapply(
    1:nrow(post),
    function(i) sum(dnorm(
      d2$height,
      mean = post$mu[i],
      sd = post$sigma[i],
      log = TRUE) )
    )

# calcul de la probabilité a posteriori (non normalisée)
post$prod <-
  post$LL +
  dnorm(x = post$mu, mean = 178, sd = 20, log = TRUE) +
  dunif(x = post$sigma, min = 0, max = 50, log = TRUE)

# on "annule" le log et on standardise par la valeur maximale (pour éviter les erreurs d'arrondi)
post$prob <- exp(post$prod - max(post$prod) )

Distribution postérieure - Grid approximation

# select random 20 rows of the dataframe 
post %>% slice_sample(n = 20, replace = FALSE)

         mu    sigma        LL      prod          prob
1  158.9950 6.939698 -1294.509 -1302.787  3.061742e-33
2  153.0653 5.180905 -1309.815 -1318.419  4.979087e-40
3  141.2060 8.899497 -1624.249 -1633.768 5.527241e-177
4  147.4372 5.834171 -1518.449 -1527.444 8.294039e-131
5  159.6985 4.703518 -1551.035 -1559.280 1.236927e-144
6  144.2211 5.005025 -2066.715 -2075.968  0.000000e+00
7  156.3819 4.226131 -1451.218 -1459.628 2.346978e-101
8  153.9698 4.778894 -1337.744 -1346.293  3.903626e-52
9  154.1709 8.396985 -1222.133 -1230.669  6.400155e-02
10 154.2714 6.638191 -1228.907 -1237.438  7.355137e-05
11 151.7588 7.919598 -1242.212 -1250.899  1.048422e-10
12 151.8593 7.140704 -1247.623 -1256.304  4.712972e-13
13 143.8191 5.557789 -1929.683 -1938.970 1.566645e-309
14 150.9548 8.849246 -1255.101 -1263.842  2.509686e-16
15 143.8191 4.527638 -2365.591 -2374.878  0.000000e+00
16 151.4573 4.577889 -1443.718 -1452.426  3.152275e-98
17 141.8090 4.201005 -3055.646 -3065.110  0.000000e+00
18 158.9950 4.552764 -1528.774 -1537.052 5.568941e-135
19 141.1055 4.326633 -3112.395 -3121.924  0.000000e+00
20 141.9095 5.783920 -2102.606 -2112.061  0.000000e+00

Distribution postérieure - Grid approximation

sample.rows <- sample(x = 1:nrow(post), size = 1e4, replace = TRUE, prob = post$prob)
sample.mu <- post$mu[sample.rows]
sample.sigma <- post$sigma[sample.rows]

Distribution postérieure - Distributions marginales

posterior_plot(samples = sample.mu, nbins = 30) + labs(x = expression(mu) )

Distribution postérieure - Distributions marginales

posterior_plot(samples = sample.sigma, nbins = 30) + labs(x = expression(sigma) )

Introduction à brms

Under the hood : Stan est un langage de programmation probabiliste écrit en C++, et qui implémente plusieurs algorithmes de MCMC : HMC, NUTS, L-BFGS…

data {
  int<lower=0> J; // number of schools 
  real y[J]; // estimated treatment effects
  real<lower=0> sigma[J]; // s.e. of effect estimates 
  }

parameters {
  real mu; 
  real<lower=0> tau;
  real eta[J];
  }

transformed parameters {
  real theta[J];
  for (j in 1:J)
    theta[j] = mu + tau * eta[j];
  }

model {
  target += normal_lpdf(eta | 0, 1);
  target += normal_lpdf(y | theta, sigma);
  }

Bayesian regression models using Stan

Le package brms (Bürkner, 2017) permet de fitter des modèles multi-niveaux (ou pas) linéaires (ou pas) bayésiens en Stan mais en utilisant la syntaxe de lme4.

model <- brm(y ~ x + (1 | subject) + (1 | item), data = d, family = gaussian() )

Rappels de syntaxe

Le package brms utilise la même syntaxe que les fonctions de base R (comme lm) ou que le package lme4.

Reaction ~ Days + (1 + Days | Subject)

mvbind(Reaction, Memory) ~ Days + (1 + Days | Subject)

mvbind(Reaction, Memory) ~ Days + (1 + Days | Subject)
mvbind(Reaction, Memory) ~ 1 + Days + (1 + Days | Subject)

Rappels de syntaxe

Si l’on veut fitter un modèle sans intercept (why not), il faut le spécifier explicitement comme ci-dessous.

mvbind(Reaction, Memory) ~ 0 + Days + (1 + Days | Subject)

brm(Reaction ~ 1 + Days + (1 + Days | Subject), family = lognormal() )

Quelques fonctions utiles

# générer le code du modèle en Stan
make_stancode(formula, ...)
stancode(fit)

# définir les priors
get_prior(formula, ...)
set_prior(prior, ...)

# récupérer les prédictions du modèle
fitted(fit, ...)
predict(fit, ...)
conditional_effects(fit, ...)

# posterior predictive checking
pp_check(fit, ...)

# comparaison de modèles
loo(fit1, fit2, ...)
bayes_factor(fit1, fit2, ...)
model_weights(fit1, fit2, ...)

# test d'hypothèse
hypothesis(fit, hypothesis, ...)

Un premier exemple

library(brms)
mod1 <- brm(height ~ 1, data = d2)

posterior_summary(mod1, pars = c("^b_", "sigma"), probs = c(0.025, 0.975) )

            Estimate Est.Error       Q2.5     Q97.5
b_Intercept 154.6045 0.4197427 153.787369 155.43024
sigma         7.7679 0.2902743   7.225066   8.36547

En utilisant notre prior

Par défaut brms utilise un prior très peu informatif centré sur la valeur moyenne de la variable mesurée. On peut donc affiner l’estimation réalisée par ce modèle en utilisant nos connaissances sur la distribution habituelle des tailles chez les humains.

get_prior(height ~ 1, data = d2)

                    prior     class coef group resp dpar nlpar lb ub tag
 student_t(3, 154.3, 8.5) Intercept                                     
     student_t(3, 0, 8.5)     sigma                             0       
  source
 default
 default

En utilisant notre prior

priors <- c(
  prior(normal(178, 20), class = Intercept),
  prior(exponential(0.01), class = sigma)
  )

mod2 <- brm(
  height ~ 1,
  prior = priors,
  family = gaussian(),
  data = d2
  )

En utilisant notre prior

summary(mod2)

 Family: gaussian 
  Links: mu = identity 
Formula: height ~ 1 
   Data: d2 (Number of observations: 352) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Regression Coefficients:
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept   154.61      0.42   153.80   155.44 1.00     3101     2416

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     7.76      0.30     7.19     8.38 1.00     3154     2598

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

En utilisant un prior plus informatif

priors <- c(
  prior(normal(178, 0.1), class = Intercept),
  prior(exponential(0.01), class = sigma)
  )

mod3 <- brm(
  height ~ 1,
  prior = priors,
  family = gaussian(),
  data = d2
  )

En utilisant un prior plus informatif

summary(mod3)

 Family: gaussian 
  Links: mu = identity 
Formula: height ~ 1 
   Data: d2 (Number of observations: 352) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Regression Coefficients:
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept   177.86      0.10   177.67   178.07 1.00     3238     2695

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    24.58      0.94    22.78    26.52 1.00     3359     2438

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

On remarque que la valeur estimée pour \(\mu\) n’a presque pas “bougée” du prior…mais on remarque également que la valeur estimée pour \(\sigma\) a largement augmentée. Nous avons dit au modèle que nous étions assez certain de notre valeur de \(\mu\), le modèle s’est ensuite “adapté”, ce qui explique la valeur de \(\sigma\)…

Précision du prior (heuristique)

Le prior peut généralement être considéré comme un posterior obtenu sur des données antérieures.

Visualiser les échantillons de la distribution postérieure

post <- as_draws_df(x = mod2) %>%
    mutate(density = get_density(b_Intercept, sigma, n = 1e2) )

ggplot(post, aes(x = b_Intercept, y = sigma, color = density) ) +
    geom_point(size = 2, alpha = 0.5, show.legend = FALSE) +
    labs(x = expression(mu), y = expression(sigma) ) +
    viridis::scale_color_viridis()